// Copyright (c) 2011 Bo Zhou<Bo.Schwarzstein@gmail.com>
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
//		http://www.apache.org/licenses/LICENSE-2.0 
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.

#include "AFMarschner.inl"

surface AFMarschner
(
	uniform float intensityD = 1;
	uniform float middleD = 0.5;
	uniform float middleRangeD = 0.1;
	uniform color colorRootD = color(1, 1, 1);
	uniform color colorTipD = color(1, 1, 1);

	uniform float intensityR = 1;
	uniform color colorR = color(1, 1, 1);
	uniform float longitudinalShiftR = -7.5;
	uniform float longitudinalWidthR = 7.5;
	
	uniform float intensityTT = 1;
	uniform color colorTT = color(1, 1, 1);
	uniform float longitudinalShiftTT = 3.75;
	uniform float longitudinalWidthTT = 3.75;
	uniform float azimuthalWidthTT = 3;
	
	uniform float intensityTRT = 1;
	uniform color colorTRT = color(1, 1, 1);
	uniform float longitudinalShiftTRT = 11.25;
	uniform float longitudinalWidthTRT = 15;

	uniform float intensityG = 1;
	uniform float azimuthalShiftG = 30;
	uniform float azimuthalWidthG = 10;
	
	uniform float attenuationFromRoot = 1.0;
)
{
	// Calcuate mixed diffuse color over curve in [0,1].
	//
	float mixingRatio = spline(v, 0, middleD, middleD + middleRangeD, 1);
	color mixedColorD = mix(colorRootD, colorTipD, mixingRatio);
	
	// Prepare vectors.
	//
	vector T = normalize(dPdv);
	vector Nn = normalize(N);
	vector B = normalize(Nn ^ T);
	vector Vn = -normalize(I);

	vector omegaO = local_spherical(Vn, B, Nn, T);

	color diffuseResult = 0;
	color scatteringResult = 0;
	illuminance(P)
	{	
		// Get \omega_o
		//
		vector Ln = normalize(L);
		vector omegaI = local_spherical(Ln, B, Nn, T);
		
		// phi is in [-PI, PI]
		//
		float phi = abs(omegaO[0] -  omegaI[0]);
		if ( phi > PI )
			phi -= 2 * PI;

		float thetaD = abs(omegaO[1] - omegaI[1]) * 0.5;
		float thetaH = (omegaO[1] + omegaI[1]) * 0.5;
		
		// Calculate the diffuse contribution.
		//
		float sinTL = sqrt(1 - pow(T.Ln, 2));
		diffuseResult += Cl * intensityD * mixedColorD * sinTL;
		
		// Get the scattering contribution.
		//
		color R = 0, TT = 0, TRT = 0;
		scatteringResult += Cl * AFMarschner(phi, thetaD, thetaH,
											 intensityR, colorR, longitudinalShiftR, longitudinalWidthR,
											 intensityTT, colorTT, longitudinalShiftTT, longitudinalWidthTT, azimuthalWidthTT,
											 intensityTRT, colorTRT, longitudinalShiftTRT, longitudinalWidthTRT,
											 intensityG, azimuthalShiftG, azimuthalWidthG,
											 R, TT, TRT);
	}

	Ci = (diffuseResult + scatteringResult) * pow(v, attenuationFromRoot);
	Oi = Os;
}

